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In this work we propose an accurate first-principle approach to calculate the transient photo¬ 
absorption spectrum measured in Pump& Probe experiments. We formulate a condition of adiabatic- 
ity and thoroughly analyze the simplihcations brought about by the fulfillment of this condition in 
the non-equilibrium Green’s function (NEGF) framework. Starting from the Kadanoff-Baym equa¬ 
tions we derive a non-equilibrium Bethe-Salpeter equation (BSE) for the response function that can 
be implemented in most of the already existing ab-initio codes. In addition, the adiabatic approx¬ 
imation is benchmarked against full NEGF simulations in simple model hamiltonians, even under 
extreme, nonadiabatic conditions where it is expected to fail. We hnd that the non-equilibrium 
BSE is very robust and captures important spectral features in a wide range of experimental con- 
hgurations. 
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I. INTRODUCTION 

The impressive progresses in ultrafast and ultrastrong 
laser-pulse technology has paved the way to the modern 
non-equilibrium (NEQ) attosecond spectroscopies. 
Unlike conventional spectroscopies, the sample is driven 
away from equilibrium by a strong laser pulse (the 
pump) before the photo-absorption of a weaker (pertur¬ 
bative) probe field is measured. Photo-absorption pump- 
and-probe (P&p) spectroscopy experiments are carried 
out using pump pulses with frequency in the infrared- 
ultraviolet range and ultrashort probe pulses (down to a 
few hundreds of attoseconds). By varying the delay be¬ 
tween the pump and probe pulses one can monitor the 
excited-state dynamics in a wide energy range. 

For samples of linear dimension (in the case of ex¬ 
tended systems this is the dimension of the primitive cell) 
smaller than the wave-length of the incident light, the 
measured signal can be calculated theoretically from the 
NEQ density response function®^^^ x (dipole approxima¬ 
tion) or, equivalently, from equilibrium dipole correlators 
of order larger than two.^®“^® In the present manuscript 
we follow the first path. 

At equilibrium x(w) can be used to construct the 
dipole-dipole correlation function a (uj) in isolated sys¬ 
tem, or the dielectric function e (uj) in extended systems. 
For correlated systems the calculation of x is, in general, 
a difficult task and one has to resort to approximations. 
The most suitable many-body scheme to implement de¬ 
pends on the sample. For atomic or small molecular sam¬ 
ples the Configuration Interaction (Cl) scheme consists in 
expanding the many-body state in Slater determinants 
to obtain eigenstates and eigenvalues. Subsequently the 
oscillator strengths are computed and used to construct 
X from a Lehmann representation. For molecules with 
tens or more nuclei as well as for crystals the number of 
Cl configurations is too large for present-day computa¬ 


tional capabilities and alternative (statistical in nature) 
approaches are required. One such approach is many- 
body perturbation theory (MBPT). In MBPT the two- 
particle electron-hole propagator L satisfies a diagram¬ 
matic equation known as the Bethe-Salpeter Equation 
(BSE) and x is constructed from a space-time contrac¬ 
tion of the arguments of The BSE has been suc¬ 

cessfully applied to study photo-absorption spectroscopy 
of systems ranging from small molecules to bulk metals 
and insulators. In this context the BSE is solved at the 
GW level with a statically screened interaction. 

Another convenient alternative to Cl (and MBPT) 
is the Linear Response (LR) Time-Dependent Density- 
Functional Theory (TDDFT).2b28 Although TDDFT 
is in principle exact,the available functionals for 
actual calculations are based on the Adiabatic Lo¬ 
cal Density Approximation (ALDA).®^^®® It is well 
known that ALDA functionals fail in capturing double- 
excitations,®^’®® charge transfer excitations®®^®® or the 
Coulomb blockade phenomenon®®’"^® in equilibrium sys¬ 
tems. For extended systems ALDA performs poorly in 
the description of the response function as it misses the 
long-range electron-hole interaction needed to describe 
excitons.®®’^® Therefore, the applicability of LR-TDDFT 
is at present restricted to weakly correlated systems with 
a spectrum dominated by single particle-hole excitations. 

Similarly to the equilibrium case the P&p photo¬ 
absorption spectrum is described by the NEQ response 
function In this work we identify a set of con¬ 

straints between characteristic times that allows us to 
rewrite x(L t') ns a function of the delay r between the 
pump and probe pulses and of the time difference t — t', 
i.e., xitjt') = X^(t ~ t')- Henceforth we will refer to 
this approximation as the adiabatic approximation. The 
mathematical rigorous definition of the adiabatic approx¬ 
imation as well as its testing in a P&p set-up is the central 
objective of the present manuscript. 
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The adiabatic response function can be computed at 
different levels of accuracy depending on the theoretical 
scheme used. In the Cl approach the time-dependent 
expansion coefficients are used to calculate the time- 
dependent product of oscillator strengths and subse¬ 
quently these products are inserted into a Lehmann-like 
representation of the NEQ adiabatic x'^(w) to yield a 
P&p spectrum with a time-dependent modulation of the 
peak intensity.®d0d2,43 MBPT, instead, we show 

that the equation of motion for x’’(< — t') can be rewrit¬ 
ten as a BSE. The main difference with the equilibrium 
BSE is that the equilibrium single-particle density ma¬ 
trix is replaced by its time-dependent value as, for in¬ 
stance, obtained from the solution of a Boltzmann-like 
equation.The NEQ adiabatic x^i^) could also be com¬ 
puted within LR-TDDET. However, it is reasonable to 
expect that the performance of ALDA functionals does 
not improve in NEQ situations. 

The structure of the paper is as follows. Section II 
presents a brief self-contained introduction to the link 
between the macroscopic observable and the microscopic 
theory. We discuss both the real-time (Section II A) 
and the response function (Section IIB) representations. 
Here we also identify a set of characteristic times in terms 
of which the condition of adiabaticity is formulated. 
The MBPT approach to x is developed in Section HI 
where we introduce the non-equilibrium Green’s func¬ 
tions (NEGP).^°’^^’^®“^® The NEGE approach is compu¬ 
tationally more expensive than TDDET but it has the ad¬ 
vantage of including dynamical correlations in a nonper- 
turbative diagrammatic fashion. To reduce the numerical 
cost we implement, in Section HI A, NEGE within the 
Generalized Kadanoff-Baym Ansatz^"'’^^ (GKBA) and 
then derive the linear response equations in Section HIB. 
Except for the GKBA no other approximations are made 
at this stage. The complexity of the problem is further 
reduced in Section IV. Here we exploit the adiabatic ap¬ 
proximation and obtain the central result of this work, 
namely a NEQ-BSE. We examine differences and analo¬ 
gies with the more standard equilibrium BSE and discuss 
the possibility of converting the NEQ-BSE into a Dyson- 
like equation in Section IV A. Einally, in Section V we 
illustrate the theory in a model system by benchmark¬ 
ing the performance of the NEQ-BSE against full NEGE 
calculations. A summary of the paper and concluding 
remarks are drawn in Section VI. 


II. THE TRANSIENT PHOTO-ABSORPTION 
SPECTRUM 

In this Section we relate the macroscopic quantity mea¬ 
sured in a P&p experiment to microscopic quantum- 
mechanical properties of the probed sample. This link 
establishes a connection between the experimental signal 
and the solution of the complex quantum kinetic equa¬ 
tion for the one-particle density-matrix. 


A. A real time approach 

In a P&p experiment the transient photo-absorption 
spectrum of a system driven out of equilibrium by a 
pump field is measured. The theoretical description of 
the driven system is achieved by evolving the many-body 
state in the simultaneous presence of the pump field and 
of a weak probe field. Let E and e be the electric pump 
and probe held respectively. We define the different terms 
constituting the many-body Hamiltonian H(t) according 
to 


ff° = f + V^, 

(la) 

= H^ + 

(lb) 


(Ic) 

H{t) = + e(t) • d. 

(Id) 


Here T is the kinetic energy operator, Vm the external 
static potential of the nuclei and 14 e the electron-electron 
interaction. Therefore is the Hamiltonian of the un¬ 
perturbed system. The inclusion of other interactions, 
e.g., the electron-phonon interaction, does not modify 
the derivation and the results of the present section. The 
terms E(t) • d and e(t) • d describe the coupling of the 
electrons with the pump and probe fields in the dipole 
approximation, d being the dipole operator (see below for 
its mathematical definition). For simplicity, we consider 
linearly polarized pump and probe fields: 

E(t) = T]pE{t), (2a) 

e{t) =’Hpeit), (2b) 

with rip and rjp the polarization vectors. The general¬ 
ization to other kind of polarizations is straightforward. 

We work in the second quantization formalism and in¬ 
troduce a suitable single-particle basis with orthonor¬ 
mal wave-functions {(pi (r)}. Then the creation and 
annihilation field-operator ■i/)^ (r) and ip (r) for a parti¬ 
cle at position r in space are expanded according to 
V’(i’) = Ei {pi(Y)ci. The one-particle density-matrix 
operator takes the form 

/5(r,r') = V>1'(r)V;(r') = Y^ipl{Y)ipj{v')pji, (3) 

with pji = cjcj . Similarly, the dipole operator projected 
along the probe field in the {pi (r)} basis reads 

d = r]p-d = J dr {'qp- r)p{r, r) = d^jpp, (4) 

with dij = f drp*(r)(qp ■ r)pj(r) the dipole matrix el¬ 
ements. In Eq. (4) and in the remainder of the paper 
we use the Einstein convention that repeated indices are 
summed over. The time-dependent expectation value of 
the dipole operator is given by 

d(t) = (4/(t)|d>(t)) = d,,(4/(t)|/5,,|vI/(Q), (5) 
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where is the state of the system at time t. 

Without any loss of generality we assume that the 
switch-on time of the pump and probe fields is larger 
than zero; hence the system is in the ground state |4'g) 
at time t = 0. Let Unit) be the unitary evolution oper¬ 
ator corresponding to a system with dynamics H 

i>0. (6) 

The time-dependent matrix elements of the one-particle 
density matrix in the presence of both the pump and 
probe fields are therefore given by 

P,,{t) = (Vl/(t)|p,,|'l-(t)) = (7) 


Replacing H with TL”®*! in Eq. (6) we have the evolution 
operator in the presence of the pump only, To 

simplify the notation we put a tilde on time-dependent 
expectation values obtained with a probe-free propaga¬ 
tion. Thus 


and hence d{t) = dijPji{t). 

For optically thin samples"'® the transmitted probe 
field e'{t) = r]pe'{t) is related to the probe-induced vari¬ 
ation (5d (t) = d(t) — d(t) of the dipole moment by'^ 

e'(t)=e(t) + ^|[<5d(t)], (9) 


where S is the cross section of the sample (assumed to 
be smaller than the cross section of the laser beam). 

The transmitted probe field is typically split in two 
halves and then merged back by a spectrometer, thus 
generating an electric field -|-e'(t — 5)] with a tun¬ 

able delay <5 > 0. In a P&p experiment the integrated 
intensity of this field, i.e., the total absorbed energy per 
unit area, is measured as a function of 6: 




pOO 

/ rif 

e'(t) -1- e'(t — (5) 

/ UjL 

/ —oo 
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( 10 ) 


The resulting function I'(5) is then cosine-transformed 


a'(w)=/ d5 I’{5 )cos{uj5), (11) 

Jo 

to gain information about the absorption energies of the 
system. Although the probe pulse has a finite duration 
the time integral in Eq. (10) goes from minus to plus 
infinity since the cosine transform requires I'{S) for all 
delays d > 0. 

Performing an analogous spectral decomposition of 
e(t) we get the intensity iJ(w) of the incident probe field. 
The photo-absorption spectrum ©(w) is therefore given 
by the difference 


©(w) 


( 12 ) 


Using Eq. (9) it is straightforward to show that'^ 

27r 

©(w) = —2oj Im [e (w) 6d (w)] — — \uj6d{uj)\^ , (13) 

be 

where e (w) and 6d (w) are the Fourier transform of 
the time-dependent probe field e{t) and probe-induced 
dipole moment 6d{t) respectively. This relation expresses 
the aforementioned link between the macroscopic inten¬ 
sity of the transmitted probe field measured in a P&p 
experiment and the microscopic dipole moment. 

In photo-absorption experiments of equilibrium sys¬ 
tems (no pump) the induced electric field [second term 
in the r.h.s. of Eq. (9)] is typically much smaller than the 
incident probe, and the quadratic term in the dipole mo¬ 
ment appearing in ©(w) can be safely discarded. More¬ 
over, Sd{uj) = a{uj)e{uj) with a{uj) = dijXji(y>)dik having 

kl 

the property that Im[a(u;)] ^ 0 for w ^ 0 (see next sec¬ 
tion), and therefore the ratio ©(u;)/|e(w)P is positive and 
independent of the shape of the probe. On the contrary, 
the photo-absorption spectrum of a pump driven system 
is not an intrinsic property of the sample since dd{Lo), 
although still linear in e, depends on e{u') at all possible 
frequencies uj'. Translating this statement from frequen¬ 
cies to times, the spectrum depends on the shape of the 
probe and on the NEQ state of the system at the time 
the probe pulse enters the sample (hence on the delay be¬ 
tween the pump and probe pulses). Furthermore, there 
might be frequencies for which the spectrum ©(w) is neg¬ 
ative due to a dominance of the stimulated emission over 
absorption. 

For a fixed shape of the pump and probe pulse the main 
interest in P&p experiments is to study the evolution 
of the spectrum as the delay r between the two pulses 
is varied. Assuming that the quadratic term in 6d, see 
Eq. (13), is small and taking into account that 6d{uj) = 
dij[pji(u!) — pji^Lo)], the resulting spectrum reads 

©"^(w) = —2a;Im [e*{u))dij5pji{u !)\, (14) 

where we define Sp = p — p. In Eq. (14) we explicitly 
added to © a dependence on t since the probe field as 
well as the time-dependent density matrix depends on 
the pump-probe delay. This dependence is, in general, 
rather complex and difficult to interpret. As we shall 
see, the calculation of the spectrum as well as its phys¬ 
ical interpretation are greatly simplified if the adiabatic 
approximation is made. 


B. A response—function representation: the 
adiabatic condition 


Equation (14) can be rewritten in a different way using 
linear response theory out of equilibrium. Let us intro¬ 
duce the (retarded) NEQ response function 


Xji it,t') = -i9 {t-t') 


Ik 


X (^sl \c\h W Cjh (i) (^') CIH it')] |4'g), (15) 


A 


U(w) — iJ^(w). 
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where Cinit) = W/jneq(t)ci Wi/noq(t) are fermion operators 
in the Heisenberg picture with respect to the probe-free 
Hamiltonian To first order in e the probe-induced 

variation of the dipole moment reads 

5d {t) = dij I dt'xji {t, t') dki e{t') 

J Ik 

= J dt'[d o ° d] e{t'). (16) 

In Eq. (16) we introduced a short-hand notation for the 
contraction of tensors of different rank. Below we define 
the four types of contractions, which include the one in 
Eq. (16), that we use in the manuscript: 

{M oV)pq = MpqVnm, (17a) 

{ToMoV) = TpqMqpVnm, (17b) 

mn 

(M O N)mn = MmnNqp, (17c) 

rs pq rs 

[N, V]mn = -[E, N]mn = (17d) 

PQ PQ pq pq 

The rank of the tensors will be clear from the context. 
Notice that Eq. (17d) has the same structure of a com¬ 
mutator since the lower indices are fixed. Taking into 
account Eq. (14) we clearly see from Equation (16) the 
relation between the P&p spectrum and the NEQ re¬ 
sponse function; we can also appreciate the complex time 
dependence introduced by the pump field. In fact, in 
equilibrium (no pump) the response function reduces to 
a function of (t — t') due to the invariance under time- 
translations. Using this invariance, the linear response 
relation Eq. (16) in Eourier space reads Sd{u!) = a{uj)e{uj) 
with a = {do xo d), and the ratio &{uj)/\e{u)\^ becomes 
independent of the probe. As already discussed in the in¬ 
troduction the equilibrium response function can be cal¬ 
culated by solving the BSE. 

In the time domain the equation for the electron-hole 
propagator L (x follows from a space-time contraction 
of L) is valid out-of-equilibrium too^*^ but its numerical 
solution is essentially impossible for present-days compu¬ 
tational capabilities. The problem is therefore to find a 
simple but still accurate approach to calculate the NEQ y 
within MBPT. Eor this purpose we will extend the equi¬ 
librium BSE to NEQ situations relevant to P&p experi¬ 
ments and provide a sound interpretation of the two-time 
dependence. In the following we refer to this equation as 
the NEQ-BSE. 

We begin the discussion by introducing two fundamen¬ 
tal characteristic times that support the adiabatic approx¬ 
imation: the key idea is that a NEQ-BSE is meaningful 
whenever the system is substantially frozen in a NEQ 
configuration during the measurement process. The char¬ 
acteristic times are 

(i) the time scale Tp of the electron dynamics induced 
by the pump. If At ^ Tp, then p{t-\- At) « p (t). 

(ii) the life-time Tp of the dressed probe pulse, which 
is the duration of the measurement process too. 


We can formulate the condition of applicability of the 
adiabatic approximation as 

Tp > Tp. (18) 

Equation (18) expresses the physical condition that the 
probe-free p{t) has to vary on a time scale (Tp) much 
longer than the duration (r^) of the dressed probe. Of 
course for Tp to be smaller than typical electronic time 
scales there should exist decay channels faster than the 
radiative decay. This is the case of solid slabs as well as of 
thick atomic or molecular gases. The following analysis 
applies to these class of systems. 

We identify two different situations where the condi¬ 
tion in Eq. (18) is fulfilled. 1) If the pump itself varies 
on a time scale Tp ^ Tp then Eq. (18) is always ful¬ 
filled since the pump-induced dynamics cannot be faster 
than Tp. In this case the adiabatic approximation, and 
hence the NEQ-BSE, can be used to describe the tran¬ 
sient spectrum for any delay t between the pump and 
probe fields. 2) In general, however, the pump is a pulse 
of duration Ap, see Fig. 1, no longer than a few hundreds 
of femtoseconds capable of inducing arbitrary fast pro¬ 
cesses. During the action of the pump the level occupa¬ 
tions change and the system polarizes. Shortly after Ap 
we have a transient period characterized by a dephasing- 
driven drop of the pump-induced polarization, we denote 
by Tpoi the polarization life-time in this nonequilibrium 
situation, and by a stabilization of the level occupations 
at some nonequilibrium value, we denote by Tcarr the 
characteristic time for the occupations to stabilize, see 
again Fig. 1. Thus, after a time Tmax = niax(Tpoi,Tcarr), 
typically Tpoi < Tcarr, we may say that the system is in a 
quasi-stationary state with carriers in some excited lev¬ 
els. In this quasi-stationary regime the time to relax 
back to the ground state is dictated by scattering pro¬ 
cesses (electron-electron, electron-phonon and electron- 
photon) and can be of the order of picoseconds. If we 
denote by Tgcatt this relaxation time-scale then we have 
Tp = Tscatt- Suppose now to probe the system in this 
quasi-stationary state with a pulse e{t) of duration Ap. 
The probe induces a polarization Sd{t) which dresses 
the bare e{t) and, in general, has a finite life-time Tp. 
Hence the duration of the dressed probe field, which co¬ 
incides with the duration of the measurement process, is 
Tp = Ap -b Tp. In this regime the condition in Eq. (18) 
is fulfilled provided that Tp is shorter than the relaxation 
time Tscatt- This is often the case as Tp is typically in the 
femtosecond range. 

In Fig. 1 we represent the dressed probe field with os¬ 
cillations of frequency Wp. Although the characteristic 
frequency Wp can be any, it is clear that it is only for 

Tp > — (19) 

Wp 

that the Fourier transform of the probe-induced dipole 
has a well defined structure in Wp. This implies that 
the life-time Tp also sets a lower limit to the frequency 
resolution of a transient spectroscopy experiment. 



5 



FIG. 1. Illustration of the characteristic times described in the main text: Tp is the time scale of the electron dynamics induced 
by the pump, Tpoi the dephasing time of the pump-induced polarization, Tcarr the stabilization time of the occupations, Tacatt 
the time to relax back to the equilibrium state and Tp the life-time of the dressed probe field. We also display the delay r 
between pump and probe. 


When the inequality of Eq. (18) is satisfied the probe 
sees a NEQ frozen system. If we take t = 0 as the time 
at which the pump is on then the probe acts at t = t and 
for times G [t — Tp^r -\- Tp] the response function 

- X'"(t - t') ( 20 ) 

depends only on {t — t') to a large extent. We will pro¬ 
vide a more precise definition of x’’)! — t') in the next 
section. For the time being we observe that whenever 
we can make the adiabatic approximation of Eq. (20) the 
transient photoabsorption spectrum of Eq. (14) can be 
written as 

&{ui) = —2u! |e(w)p Im [d o x’^(a;) o d]. (21) 

Consequently the ratio 6'^(a;)/|e(a;)p becomes indepen¬ 
dent of the probe and can be interpreted as an intrinsic 
property of the nonequilibrium system. 

As a very general remark we notice that when the sys¬ 
tem is probed after the pump (no overlap between the 
pulses) the probe-induced dipole moment oscillates at 
frequencies Xla /3 = Ea — Ejs, where Ea,Ep are eigen- 
energies of see Eqs. (15) and (16).^^’®^ Furthermore, 
the amplitude of the oscillations depends on the delay 
t}'^ Therefore P&p spectra are richer than equilibrium 
spectra where the probe-induced dipole moment can os¬ 
cillate only at frequencies Xlp = Eg — Ep, with Eg the 
ground state energy, with constant amplitudes. The ex¬ 
tra transitions are usually referred to as photo-induced 
absorption and stimulated emission. 

To summarize, Eqs. (14) and (21) represent two dif¬ 
ferent ways of calculating the transient photo-absorption 
spectrum. We could either perform a time propaga¬ 
tion with both the pump and the probe, a second time 
propagation with only the pump and then extract the 
probe-induced dipole moment 6d or we can evaluate the 
response function x^ from a NEQ-BSE. The latter ap¬ 
proach is developed in the next section. 


III. A NON-EQUILIBRIUM GREEN’S 
FUNCTION APPROACH TO TRANSIENT 
ABSORPTION 

In the previous section we have introduced the theoret¬ 
ical description of transient absorption experiments with 
two possible approaches. The first, which is exact, based 
on Eq. (14) and the second, which uses the adiabatic 
approximation, based on Eq. (21). However, these equa¬ 
tions assume that it is ideally possible to compute the 
exact time-dependent density matrix or the exact adia¬ 
batic response function. This is not doable in practice 
and one has to resort to approximations. In the follow¬ 
ing we show how to use NEGF theory to obtain a MBPT 
equation for 5d (t ). In the next Section we use this result 
to generate an equation for the NEQ response function 
x(t, t'), and subsequently make the adiabatic approxima¬ 
tion to derive the NEQ-BSE for x^(i — t'). 

In the MBPT approach the description in terms of the 
many-body hamiltonian containing the electron-electron 
interaction is replaced by a description in terms of the one 
particle hamiltonian and the many-body self-energy. We 
thus define: 





VNir) 


= h° + S°, 

h*(t) = -t AS* -t E(<) • d, 
h*(t) = h*{t) -\- (5E* -I- e(f) • d. 


(22a) 

(22b) 

(22c) 

(22d) 


We have here introduced the t superscript to indicate a 
quantity whose time-dependence is given by the implicit 
dependence on the density matrix. This means that, for 
a generic function /, we have 


(23) 

f{t)^f[p{t)]{t). (24) 
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The difference between /* and /* (t) is that the second 
function has an explicit time dependence too. Further¬ 
more, to indicate that the function is calculated at the 
probe-free density matrix p we put a tilde symbol on 
the function. Thus f* = f[p (t)] and /*(t) = f [p (t)] (t). 
Let us define the three different self-energies appearing 
in Eqs. (22). The self-energy is the static 

part of the equilibrium many-body self-energy and it 
is therefore calculated at the equilibrium density matrix 
p®*!. The self-energies AS* and JS* are the variations 
due to a change in p induced by the pump and the probe 
respectively: 

AS* = S,[p(t)] - S,[p®'i] = S* - SO, (25a) 
5S‘ ^ S, [pit)] - S, m] = - S*. (25b) 

In general S^ is the Hartree-Fock (HF) plus static cor¬ 
relation self-energy. It plays a crucial role as it renormal¬ 
izes the single-particle level energies and introduces cor¬ 
relation effects (like electron-hole attraction) also in the 
polarization function. The different possible approxima¬ 
tions to S^ reflect the different kind of physics introduced 
in the dynamics: 

(i) A mean-field potential that mimics the correlation 
effects. An example is DFT where S* is local in 
space and given by the sum of the Hartree and 
exchange-correlation potential. 


It is easy to verify that the one-particle density matrix 
is given by the lesser Green’s function at equal times, 
p it) = —iG^ it, t). The functions satisfy a set of cou¬ 
pled equations known as the Kadanoff-Baym equations 
The KBE are integro-differential 
equations with a self-energy kernel depending on both 
G< and G>. It is possible to collapse the KBE into 
a single equation for the one-particle density matrix by 
making the so called Generalized Kadanoff-Baym Ansatz 
(GKBA).^° The corresponding equation for p reads 

±pit) + ^[h\t),pit)]=-I\t), (28) 

where Ii*(t) is defined in Eq. (22). 

The collision integral /*(t) = /[p(t)](t) on the r.h.s. 
of Eq. (28), at difference with the static self-energies 
previously discussed, is non-local in time (unless spe¬ 
cific approximations are made). The functional form is 
uniquely determined through the GKBA once an approx¬ 
imation for the correlation self-energy, Ec, is made. Let 
us show how to obtain /* starting from its exact KBE 
expression and then making the GKBA. From the KBE 
we have 


lit) = / dt 


^<it,t)G^-\i,t) 

+EW(t,t)G<(t,t) 


-f H.c. (29) 


(ii) HF self-energy. In this case no correlation is in¬ 
cluded. The HF self-energy reads E* = K o pit), 
with the four-index tensor Vij = 2vimnj — Vimjn 

mn 

and Vimnj the two-electron Goulomb integrals. 

(iii) Hartree plus a Goulomb Hole and Screened Ex¬ 

change (GOHSEX) self-energy. In this case corre¬ 
lation is included using a linear-response approx¬ 
imation but dynamical effects are neglected. The 
GOHSEX self-energy reads E* = K* o p -j- with 
VG = 2vimnj - ^^d the Goulomb hole 

mn 

potential. In H* the screened exchange interaction 
reads 

v*(r,r') = J dreRpA[p(t)](r,r)u(r-r'). (26) 

In all cases the static self-energy is a time-local func¬ 

tional of the density matrix. 

A. Real-time dynamics (I): the Generalized 
Kadanoff-Baym Ansatz 

In NEGF theory the key quantities are the lesser, 
G< it, t'), and greater, G^ (t, t'). Green’s functions. These 
functions are defined according to 

= f(c]^(t')ciff(0). (27a) 

= -*(ciff(f)c]ff(t'))- (27b) 


with Ec a functional of G^ and G^. The functional 
form of Ec must be consistent with the choice of Eg, 
i.e. Ec = E — Eg with E the full many-body self-energy. 
Retarded/advanced functions carry a superscript (r)/(a) 
and are defined in terms of the lesser and greater func¬ 
tions according to 

AW(t,t') = [xG)it',t)]^ 

= eit -1') [x> it, t') - x< it, t')], (30) 

where X can be G, Ec or any other two-time correlator. 
The GKBA is an ansatz for G^ which turns Ec, and 
hence the collision integral, into a functional of p and 

GW/(a). 

G<it,t') = - G^'^Ht,t')pit') + pit)GG')it,t'), (3Ia) 
G>it,t') = + G^'^Ht,t')pit')-pit)GGlit,G), (31b) 

where p = 1 — p. To transform /(t) into a functional 
I* it) of the density matrix, and hence to close Eq. (28), 
one needs to express the propagator in terms of p. 
Depending on the system there exist optimal approxima¬ 
tions to the propagator, the most common one being the 
quasi-particle (QP) propagator 

GW it, t') = -ieit - t')Te-^ . (32) 

For (small) finite systems the choice (usu¬ 

ally is the HF single-particle hamiltonian) is a good 
choice. For extended systems, however, the lack of damp¬ 
ing in /i®'^ prevents the system to relax. In these cases 
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the propagator is typically corrected by adding non- 
hermitean terms given by the quasi-particle life-times 
h^P = + ^ 7 . 54-59 


B. Real-time dynamics (II): the linear regime 

If the probe is a weak perturbation we can work within 
a linear-response approach. Then (5E* is of first order in 
e and the collision integral can be expanded as 

l\t)^P{t) + 6P{t). (33) 

Inserting Eq. (33) into Eq. (28) and equating terms of 
the same order in the probe field we get two equations, 
one for p and another for 5p (omitting the explicit time- 
dependence from the various quantities): 

—p + i\h*,p]=—P, (34a) 

— (5p -I- i[h^ ,5p] + [(5E* -f e • d,p] = —51*. (34b) 

As we are in the linear-response regime, we can rewrite 
(5S‘ and 51* in terms of kernel functions of the probe-free 
density matrix p. The notation introduced proves now 
useful because it highlights the dependence on p and 5p-. 

5^l=Klo5p{t), (35a) 

51* (t) = ( diK*{t,t) o 5p{t). (35b) 


The static kernel AT* and the correlation kernel K* de¬ 
pend only on p. Furthermore Kl(t,t) vanishes for i > t 
since I*{t) depends on p(t) only for t < t, as it follows 
directly from Eq. (29) and the GKBA in Eqs. (31). With 
Eqs. (35) we can rewrite Eq. (34b) as 

^5p + i[h*,5p] + i[Klo 5p + e ■ d,p] = 

— j dt o 5p{t). (36) 


Equation (36) is the many-body equation for the calcu¬ 
lation of the probe-induced change of the density matrix. 
In the next section we combine Eq. (36) with the condi¬ 
tion of adiabaticity in Eq. (18) to derive a NEQ-BSE. 


IV. NON-EQUILIBRIUM BETHE-SALPETER 
EQUATION 

The next step in the derivation of a BSE in the presence 
of the pump field is to transform Eq. (36) into an equation 
for the response function. To this end we use the relation: 


with Uki(t) = e{t) ■ dki- Taking the functional derivative 
of Eq. (36) with respect to u{t') and find 




K*s o + 15it - t'), p{t) 


- / dt K*pt,t)ox{t,t'), (38) 


where we introduced the four index tensor \ji = 5ji5ik. 

kl . 

At zero pump this equation reduces to the equilibrium 
BSE 

+ ^ [K!^ o x^pt - t') + 15it - t'), ^*=9] = 

- J dt (39) 


The differences between Eq. (38) and Eq. (39) are: 

(i) in the equilibrium limit all quantities depend on the 
relative time coordinate only; 

(ii) is time independent while p{t) is time depen¬ 
dent; 

(iii) the static equilibrium hamiltonian /i®9 is replaced 
by the time dependent h*{t)] 

(iv) the kernels A* and K*{t,t) are evaluated at 

the pump-driven time-dependent density matrix p 
whereas the kernels Ar|9 and — t) are evalu¬ 

ated at the static equilibrium density matrix /j®®*. 

Due to these points it is not possible to reduce Eq. (38) to 
an algebraic equation for x(t,t'), as it is commonly done 
in state-of-the-art equilibrium calculations after Fourier 
transforming with respect to the time difference t — t'. 
In Eq. (38) x(t, t') is not a function oit — t' and, further¬ 
more, the dependence on t appears both implicitly and 
explicitly in h*, Kl, K* and p. 

Analytical progress can be made provided that the adi¬ 
abatic condition, see Eq. (18), is fulfilled. We recall that 
in this approximation the pump-driven density matrix 
p{t) varies slowly over the life-time Tp of the dressed 
probe. Thus for t G [t — Tp,T + Tp] we have 

p{t)^p{T). (40) 

In the same time window E(<) « E(t) and hence Eq. (40) 
implies that 


h*it) « h^{T), (41) 

Kl « Kl. (42) 


Ik 


Sp^t) 

5uki{t'Y 


( 37 ) 


This is a direct consequence of the fact that the function¬ 
als h*{t) and AT* are time-local functionals of p(t). 



Another simplification brought about by the adiabatic 
condition is that for times t,t' G [t — Tp,T + Tp] the re¬ 
tarded Green’s function, see Eq. (32), can be approxi¬ 
mated as 

ss —— t') exp[—— t')]. (43) 

Therefore, the adiabatic retarded Green’s function is in¬ 
variant under time translations. The crucial consequence 
of this fact is that the correlation kernel too becomes a 
function of the time difference only: 

(44) 

Taking into account Eqs. (40-44) we see that the solu¬ 
tion of Eq. (38) is a response function x(t, t') ~ — t') 

depending on the delay r and on the time difference t — t'. 

In the adiabatic approximation Eq. (38) can be con¬ 
veniently Fourier transformed to yield an algebraic equa¬ 
tion for the frequency-dependent response function 

-iujx''iuj)+i h'"(T),x^(a;) +i kj o x^( uj) + 1, pir) 

= -k^{uj)ox'"{uj). (45) 

This is the aforementioned NEQ-BSE and the main re¬ 
sult of the present work. We emphasize that x’’ is the 
response function of the hnite system. In the case of 
extended systems x’" is equivalent to the macroscopic 
response function obtained from a supercell calculation 
where the spatial long-range component of the induced 
Hartree held (corresponding to its q —)■ 0 Fourier compo¬ 
nent) has been removed^^. 

The solution of Eq. (45) requires a preliminary calcu¬ 
lation of the one-particle density matrix p{t). In the next 
sub-section we show how to rewrite the NEQ-BSE as a 
Dyson equation for x^ ■ The NEQ Dyson equation is then 
compared with its equilibrium counterpart to provide an 
intuitive physical interpretation of the response function. 

A. Reduction to a Dyson equation 

The NEQ-BSE, Eq. (45), can be implemented in most 
of the ab-initio numerical schemes and codes. However, 
in order to create an even closer connection to standard 
implementations of the BSE we further discuss the ap¬ 
proximations and conditions under which Eq. (45) turns 
into a simple Dyson equation. 

The crucial aspect is the choice of the reference ba¬ 
sis and its link with the adiabatic approximation. Let 
us first re-examine the equilibrium case. Gonsider the 
representation in which is diagonal, i.e., 

Then, the equilibrium density matrix is diagonal too and 
its entries are the occupation factors of the electronic lev¬ 
els: p®^ = Sijk^. In this basis the Fourier transform of 
the equilibrium BSE, i.e. Eq. (39), reads 

[wl - Ae®‘i + iK^^iio)] o x®'^(w) = 

-Aroil+Kt^ox^^k^)], (46) 


where 

(Ac®'!),, =(er-e;'^)l.i, (47) 

pq pq 

and 

(48) 

pq pq 

Introducing the response function 

Xo‘^(w) = - [wl - Ae®'i + iK^^iuj)]-^ o Af^ (49) 

we can rewrite Eq. (46) in the form normally used in first 
principles calculations 

X^kco) = xlk^) + X?i^) o o x^kco) . (50) 

The correlation kernel appearing in Xo*^ deserves 

a comment. In most of the applications is usu¬ 
ally replaced by a constant, i.e. ss p. More 

sophisticated approximations with (AT®'^) ij ss pkij + 

mn 

'Im.Smn, have been explored.®"^ In this case the quasi¬ 
particle line-widths "ft are calculated from equilibrium 
MBPT. The approximation of a static correlation kernel 
is based on the observation that dynamical corrections 
to the screened interaction are partially cancelled by the 
dynamical effects in the quasi-particle corrections, see 
Ref. 55. 

Let us now consider the NEQ-BSE, i.e. Eq. (45). Like 
in the equilibrium case we would like to introduce a xo 
and turn Eq. (45) into a Dyson equation. However, in 
the NEQ case neither p(r) nor h'’’(r) are diagonal in the 
eigenbasis of Of course we can rotate the equilib¬ 
rium basis so to have k^r) diagonal but, in general, p ( t ) 
has off-diagonal entries in this new basis too. Let 0{t) 
be the orthogonal matrix of the transformation from the 
equilibrium basis to the adiabatic basis in which (t) is 
diagonal 

0\T)K^0{r) =6ijei{T). (51) 

L J ij 

The NEQ-BSE Eq. (45) in the adiabatic basis reads 

wl — Ae(r) -I- iAJ(w) ox’^(w) = 

- pij),klox'k) + '^ ( 52 ) 

where the four-index tensor Aeir) is defined as in 
Eq. (47) with e®'^ —>• ei(T). Next we dehne the NEQ 
response function Xo according to 

X5(w) = - w - A7(r)-kiAc (‘^) o[p('r),l] (53) 

which generalizes Eq. (49) to nondiagonal density matri¬ 
ces. Using the identity 

ox'^iuj) =[p{T),l]okl ox^k) (54) 
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we can rewrite the NEQ-BSE in a Dyson-like form 


V. A NUMERICAL EXAMPLE 


x^(w) = xl{^) + xl{^) ° Kl ox'"(w). (55) 

The analogy between Eq. (55) and the standard, equi¬ 
librium BSE becomes more evident if we make some fur¬ 
ther approximations that are often used in actual imple¬ 
mentations: 

(i) constructed from the dynamical GW self¬ 
energy. In this way the equilibrium basis is the 
quasi-particles basis whose states are renormalized 
by dynamical effects. The same approximation ap¬ 
plies, for internal consistency, to ASJ. 

(ii) In order to recover the equilibrium limit of standard 
BSE implementations the term must be the 
statically screened COHSEX approximation. This 
has been proved in Ref. 60. 

With these two approximations in mind we discuss 
Eq. (45) in the case of a weak pump Held. This condi¬ 
tion is often realized in P&p experiments as it allows to 
photo-excite the system without changing too much its 
electronic and optical properties. Therefore, weak pump 
fields provide a non-invasive method to monitor the ex¬ 
cited states of the equilibrium system. 

For weak pump fields the density of the excited carri¬ 
ers is small. This implies that we can approximate the 
orthogonal matrix Oijir) « Sij. In other words the adi¬ 
abatic basis and the equilibrium basis are essentially the 
same. The obvious and physically intuitive consequence 
of this fact is that the diagonal elements of the density 
matrix are the NEQ occupations pair) = fiij) whereas 
the off-diagonal elements describe the polarization of the 
system. If the photo-excited carrier density is small the 
off-diagonal elements can be neglected and Eq. (53) sim¬ 
plifies to 


A. Model 

We illustrate the theory developed in the previous sec¬ 
tions by calculating the transient photoabsorption spec¬ 
trum of a four-level model system with two valence states 
(orbital quantum numbers p = 1,2) and two “conduc¬ 
tion” or excited states (orbital quantum numbers p = 
3,4). In second quantization the equilibrium Hamilto¬ 
nian reads 

H ^ = 'y ( f-fiflfia + 2 'y ' (58) 

fl(T 

(7<t’ 

with the occupation operator of level p 

with spin a. The system is driven out of equilibrium by 
a strong pulse which pumps electrons from the valence 
states to the conduction states. In accordance with the 
notation of Eqs. (1) we consider a pump-dipole coupling 
of the form 

E(t) • d = E{t) ^ ^ + H.c.), (59) 

— <7 

iv=3,4 

where = {Lp^\rip-v\{p^) = ryp-d^^. After a time r the 
excited system is irradiated by a weak ultrafast probe. 
The probe-dipole coupling is the same as in Eq. (59) 
except that the field amplitude E{t) is replaced by the 
amplitude e{t) of the probe pulse. For the numerical 
simulations we choose the amplitudes E{t) and e{t) as^° 

E(t) = Eq sin^ (e~^) sinwpt (60) 

Ap 

for 0 < t < Ap and zero otherwise, and 


x5(w) 


w — Ae(r) -I- iK^ (w) 


o A/(t) 


(56) 


e(t) = Co sin^ (tt --) sina;p(t — Ap — r) (61) 

ZAr) 


where the four-index tensor A/(r) is defined as in 
Eq. (48) with —)• /^(t). In addition, the pump- 
induced renormalization of the single-particle energy lev¬ 
els is 


for0<t — Ap — T< Ap and zero otherwise. 

The equation of motion for the single-particle density 
matrix is Eq. (28). For h*{t) we take the HE Hamiltonian 
(see discussion just before Section HI A) 


G(r)«C + E(r)-d,, + AEW, (57) 

where AEJ is a time-local functional of the occupations 
only. The r-dependent renormalization of the energy lev¬ 
els represents the explanation in MBPT language of the 
well known band gap renormalization effect, i.e., the re¬ 
duction of the elemental gap induced by pump excited 
carriers. 

Another consequence of the diagonal structure of the 
density matrix is that the static kernel too becomes a 
functional of the NEQ occupations only: Kl « Ks[f{T)]. 
This dependence can be used to interpret the renormal¬ 
ization of the electron-hole interaction and hence, in sys¬ 
tems with bound excitons, the renormalization of the ex- 
citonic binding energy. 


[c^ T y ^ 2^paPQa(Q] pi/ 

Ol 

+ [E{t) + e{t)]dp^. (62) 

For the collision integral we consider a two-step relax¬ 
ation approximation (in matrix form) 

i\t)^{TP°\t),pit)-p^^} + { it), p{t)-p‘^^}, 

(63) 

where the curly brackets signify an anticommutator. In 
Eq. (63) the first term accounts for the dephasing of the 
pump-induced polarization and is responsible for driving 
the system toward a quasi-stationary state described by 
After the dephasing rP°^(t) « 0 and the collision in¬ 
tegral is dominated by the second term which describes 
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(a) 


(b) 



FIG. 2. Upper panel: (a) non-interacting energy levels ei = 0, 
£2 = 0.1, £3 = 1.0, £4 = 1.3 and (b) HF energy levels 


- Sfj,aVfii_,) with = 0.4, Wl2 = W23 = 0.2 
and otherwise = 0.1 (we recall that = p^a)- The 
position of the poles of the equilibrium TD-HF 
(solution of Eqs. (49-50) with = 0) can be calculated 
analytically and are also indicated. Lower panel: Plot of the 
damping functions 7 i°*(t) (reduced by a factor of 

10) and 72 °*(f) used in the numerical simulations. Energies 
£^, £™ and are in eV while 7 ®'^'^** and 75 '_ 2 ^ are in meV. 



FIG. 3. From bottom to top, pump pulse E{t) (black), probe- 
free dipole d{t) (red), occupation ns — n 4 = n oi valence 
states (blue) and probe-induced dipole 5d{t) (green) for the 
large 7 ^°* (top panel) and small 71 °* (bottom panel) damping 
functions. These results are obtained for a delay r = 130 fs. 
The quantities E{t), d{t) and 5d{t) are in arbitrary units. 


the relaxation toward the equilibrium state. The damp¬ 
ing matrices rP°'(t) = 7 P°*(t)l and r"°’^“(t) = 7 "““(t)l 
are proportional to the identity matrix, thus guaran¬ 
teeing the conservation of the total number of particles 
N = 2Tr[p]. Since there is no pump-induced dephasing 
in the absence of the pump rP°* is proportional to the 
amplitude of the pump pulse. 

The system has filled valence states and empty con¬ 
duction states at time t = 0, hence = diagjl, 1 , 0 , 0 }. 
The model parameters as well as the HF equilibrium con¬ 
figuration can be found in Fig. 2. For the dipole matrix 
we use 


/ 0 0 1 1 \ 
I 0 0 1 1 I 
110 0 
\ 1 1 0 0/ 


(64) 


The damping functions and the two different yP®* = 
7 ^°*, 72 °' that we consider are illustrated in the lower 
panel of Fig. 2. In particular yP®* is responsible for 
the relaxation toward the quasi-stationary density ma¬ 
trix = diag{0.9, 0.9,0.1, 0.1}. For the external fields 
we study a pump pulse of duration Ap = 66 fs and fre¬ 
quency ujp = 0.6 eV, and a probe of duration Ap = 20 


fs and frequency Wp = 0.6 eV; the amplitudes Eq, cq and 
do are chosen to yield Eodo = 0.1 eV and eodo = 0.001 
eV. 

We calculate p{t) in the presence of both pump and 
probe as well as the probe-free p and then extract 
the probe-induced dipole moment Sd{t) = d^^[p,j^{t) — 
Pv 4 i{t))- Successively, we obtain the transient spec¬ 
trum of Eq. (14) by Fourier transforming the function 
5d{t) X with fp = 80 fs the life-time of the probe- 

induced dipole. In the figures below the exponential 
damping is always included in the probe-induced dipole. 
The probe-free p is also used in Eq. (45) to calculate the 
adiabatic NEQ response function and hence the transient 
spectrum according to Eq. (21). The quality of the adi¬ 
abatic approximation is assessed in different regimes. 


B. Results and discussion 

In the upper panel of Fig. 3 we show some relevant 
quantities obtained from the numerical solution of Eq. 
(28) with yP®' = 7 }®', namely (from bottom to top) the 
pump pulse E{t), the probe-free dipole d{t), the time- 
dependent occupation no{t) = n^lt) = n(t) of the va- 
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lence states 3 and 4, and the probe-induced dipole 6d{t). 
The behavior of these quantities resemble the behavior in 
Fig. 1. When the probe arrives (t = 130 fs), the pump- 
induced polarization is completely dephased (rpoi ~ 100 
fs), and the system is slowly moving around the quasi¬ 
stationary excited state described by In this situ¬ 
ation the time-scale over which the one-particle density 
matrix changes is Tp ^ x/^scatt ^ fs. Since Tp is 
much larger than the life-time Tp = fp + Ap = 100 fs of 
the dressed probe the adiabatic condition is fulfilled, see 
Eq. (18). 

The transient absorption spectra ©’’(w) obtained 
within NEGF according to Eq. (14) and with the NEQ- 
BSE according to Eq. (21) are displayed in Fig. 4. As ex¬ 
pected, the NEQ-BSE approach is very accurate for de- 
lays T < —Tp and r > Tpoi, i.e., when the probe-induced 
dipole does not overlap the pump-induced polarization. 
For T < —Tp the spectrum 6'^(oj) is the equilibrium spec¬ 
trum with four peaks at energies uj^^, = + Aw 

(with /i = 1 , 2 , ly = 3,4 and Aw = 3 ui 3 — 2vi2 — vii), 
thus NEGF and NEQ-BSE obviously agree. For r > Tpoi 
the system is in a nonequilibrium state and the condition 
of adiabaticity matters. The NEQ-BSE well captures the 
T-dependent structure of the NEGF spectrum, with the 
correct bending of the position of the four main absorp¬ 
tion peaks towards their equilibrium value for large r. 
At first sight, the agreement seems rather good in the 
overlap region —Tp < t < Tpoi too. However, in this re¬ 
gion ©'^(w) is very small due to the sizable broadening 
induced by the large 7 f°*, and a more careful comparison 
between the NEQ-BSE and NEGF spectra reveals some 
discrepancies (not shown). 

A second simulation has been performed using the 
smaller damping function 7 P°* = 72 °^- In this case the 
pump-induced polarization d(t) is long-lived, as it can be 
seen in the lower panel of Fig. 3. After a time r = 130 
fs the time-scale Tp over which the one-particle density 
matrix changes is given by the period of the coherent os¬ 
cillations of d{t) and it is roughly equal to the inverse 
gap l/w 23 ~ 10 fs. Thus the condition of adiabatic¬ 
ity Tp ^ Tp is not fulfilled and no agreement between 
NEQ-BSE and NEGF is expected. The transient ab¬ 
sorption spectra are displayed in Fig. 5 showing that the 
two approaches differ whenever the probe experiences a 
sizable d{t), i.e., for —Tp < t < 500. In this region the 
NEGF spectrum exhibits alternating fringes character¬ 
ized by a large oscillation of the spectral weight at fixed 
w as a function of t. These features origin from the 
nonadiabatic coherent motion of the electrons between 
valence and conduction states, and hence they are out of 
reach of the NEQ-BSE approach. Remarkably, however, 
the NEQ-BSE captures important spectral features even 
in this strongly nonadiabatic situation, the most promi¬ 
nent feature being the upward bending of the main peaks 
around t = 0. When the coherence is destroyed by the 
dephasing, i.e., for t > 500 fs, the NEQ-BSE and NEGF 
spectra are found to be in excellent agreement. 


0. 0.2 0.4 0.6 0.8 1 



NEGF yi””’ 



- 150-100 -50 0 50 100 200 500 1000 

T (fs) 

0. 0.2 0.4 0.6 0.8 1 __.^pol 


NEQ BSE y; 



T (fs) 


FIG. 4. Transient absorption spectrum S’^(w) (normalized to 
its maximum value) obtained within NEGF according to Eq. 
(14) (upper panel) and with the NEQ-BSE according to Eq. 
( 21 ) (lower panel) using the damping function 7 ®’°* = 7 i°*. 
The peaks of the NEQ-BSE spectrum have been broadened 
by the inverse life-time 1/fp of 5d{t). Notice that the scale 
of the horizontal axis is linear for t < 100 fs and logarithmic 
otherwise. 


VI. CONCLUSIONS 


We propose a practical method based on MBPT to cal¬ 
culate P&p spectra for delays in the “adiabatic” regime. 
Starting from the KBE for the Keldysh Green’s function 
we use the GKBA to obtain an equation of motion for 
the one-particle density matrix p in the presence of both 
pump and probe fields. Linearization around zero probe 
yields an equation for the NEQ response function y(f, t'). 
After the action of the pump we identify a physically rel¬ 
evant regime during which the probe-free density matrix 
p varies on a time-scale much longer than the life-time of 
the dressed probe. In this regime we make the adiabatic 
approximation and show that y(f, f') can be written as 
a function of the pump-probe delay t and of the relative 
time {t — t'), i.e., x(t,t') « x'”(t — t'). This simplihca- 
tion allows us to Fourier transform with respect to the 
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FIG. 5. Same as Fig. 4 except that the small damping func¬ 
tion 7 |’°* has been used in Eq. (63). 


relative time and to derive the main result of this work: 
a NEQ-BSE that can be implemented in most of the 
ah-initio numerical schemes and codes. We further pro¬ 
vide a sound physical interpretation of the NEQ response 
function and showed that it can be related to intrinsic 
spectral properties of the nonequilibrium system. Well 
known effects like the renormalization of the band-gap 
and excitonic binding energies in semiconductors and in¬ 
sulators are naturally explained. 


The computational advantage of the NEQ-BSE over 
NEGF simulations is enormous as only the probe-free 
one-particle density matrix p enters in the solution of the 
NEQ-BSE. This implies that a single time-propagation 
is sufficient to obtain the transient spectrum for several 
delays. In contrast, the NEGF approach requires a time- 
propagation for every delay (to obtain the one-particle 
density matrix with pump and probe fields) in addition 
to the time-propagation to obtain p. The validity of 
the NEQ-BSE has been successfully demonstrated in a 
simple four-level model system and it is currently under 
investigation in more realistic hamiltonians with encour¬ 
aging results. 
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